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ABSTRACT 

We present results from numerical simulations of the interaction of internal gravity 
waves (IGW) with magnetic fields in the radiative interior of the Sun. In this second 
paper, the waves are forced self-consistent ly by an overlying convection zone and 
a toroidal magnetic field is imposed in the stably stratified layer just underneath 
convection zone. Consistent with the results of previous analytic and simple numerical 
calculations, we find a strong wave-field interaction, in which waves are refiected in the 
field region. The wave-field interaction and wave refiection depend on the field strength 
as well as adopted values of the diffusivities. In some cases wave reflection leads to 
an increased mean flow in the fleld region. In addition to reproducing some of the 
features of our simpler models, we find additional complex behaviours in these more 
complete and realistic calculations. First, we find that the spectrum of wave generation, 
both in magnetized and un-magnetized models, is not generally well described by 
available analytic models, although some overlap does exist. Similarly, we find that 
the dissipation of waves is only partially described by the results of linear theory. We 
find that the distortion of the field by waves and convective overshoot leads to rapid 
decay and entrainment of the magnetic field which subsequently changes the wave-field 
interaction. In addition, the field alters the amount of wave energy propagating into 
the deep radiative interior, at times increasing the wave energy there and at others 
decreasing it. Because of the complexity of the problem and because the durations of 
these simulations are shorter than the anticipated timescale for dynamical adjustment 
of the deep solar interior, we are unable to draw a definitive conclusion regarding the 
efliciency of angular momentum transport in the deep radiative interior by IGW in 
the presence of a magnetic field. 



1 INTRODUCTION 

It has long been known that internal gravity waves (IGW) 
can transport angular momentum over long distances. In 
the Earth's equatorial stratosphere the dissipation of IGW 
are the primary driving mechanism of the Quasi-Biennial 
Oscillation (QBO) (Baldwin et al. 2001) in which the mean 
zonal flows oscillate with a period of approximately 27-28 
months. The ability of IGW to drive mean flow oscillations 
and hence, its relevance to the QBO, was demonstrated in 
the experiment by Plumb & McEwan (1978). 

Because of their ability to transport angular momentum 
over long distances and their anticipated presence in stellar 
radiative interiors, there is considerable interest in IGW in- 
duced transport and mixing processes in stellar interiors . 
Nearly 30 years ago Press (1981) studied the propagation 
and dissipation of IGW in radiative interiors, showing that 
these waves could possibly lead to increased opacity and 
species mixing. Garcia-Lopez & Spruit (1991) (hereafter re- 
ferred to as GLS) showed that mixing by IGW could con- 
tribute to the enhanced Lithium depletion observe in main 
sequence F-type stars. Schatzman has contributed numer- 



ous studies on the effects of IGW in stellar interiors. He has 
shown that their properties allow a host of transport possi- 
bilities, such as the ability to induce diffusion and hence mix 
species (Schatzman 1996; Montalban & Schatzman 2000, 
1996), transport angular momentum (Schatzman 1993) and 
even affect opacities and hence, the solar neutrino problem 
(Schatzman 1999). 

More recently, internal waves in stellar radiative inte- 
riors have received signiflcant attention in connection with 
claims (Kumar & Quataert 1997; Zahn et al. 1997) that they 
might be the mechanism responsible for the uniform rota- 
tion of the solar interior, as inferred by helioseismic inversion 
(Thompson et al. 1996, 2003). These claims recieved heavy 
scrutiny (Ringot 1998; Gough & Mclntyre 1998) because of 
the known propensity for IGW to force shear flows rather 
than uniform flows (Plumb k, McEwan 1978; Baldwin et al. 
2001). More recent incarnations of the IGW theory for the 
uniform rotation of the solar radiative interior have incorpo- 
rated an oscillating shear layer, coined the SLO (shear-layer 
oscillation) as a filter on IGW propagating to the deep in- 
terior (Kumar et al. 1999; Talon et al. 2002; Charbonnel & 
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Talon 2005) (we will hereafter refer to Kumar et al. (1999) 
as KTZ). However, the validity of these models has also been 
questioned (Rogers & Glatzmaier 2006; Rogers et al. 2008; 
Denissenkov et al. 2008), again because the general propen- 
sity of IGW to enhance shear layers. In all of these works the 
main uncertainty arises from our incomplete understanding 
of the spectrum and amplitudes of IGW generated. 

While there have been numerous studies on the prop- 
agation, dissipation and transport properties of IGW in 
stellar interiors, very few (Schatzman 1993; MacGregor & 
Rogers 2010; Rogers & MacGregor 2010) have addressed 
the issue of magnetism. The main mechanism thought to 
generate IGW in the solar radiative interior is convective 
overshoot at the base of the convection zone. This is the 
likely storage region for the strong toroidal magnetic field 
which gives rise to the sunspot cycle. Indeed, the sites of 
IGW generation and magnetic field storage are likely to co- 
incide in many stars. As pointed out by Schatzman (1993) 
a magnetic field could prevent the propagation of IGW into 
stellar interiors. For these reasons it is crucial to understand 
the generation, propagation and dissipation of IGW in the 
presence of magnetic field. We began this study with ana- 
lytical models (MacGregor & Rogers 2010) and simple nu- 
merical simulations (Rogers & MacGregor 2010). Here we 
extend those studies to more realistic numerical simulations 
in which the IGW are self-consistently driven by an overly- 
ing convection zone. We then study the interaction of these 
waves with an imposed toroidal field. In Section 2 we de- 
scribe our numerical model, in Section 3 we discuss the gen- 
eration and propagation of IGW in the presence (3.2) and 
absence (3.1) of a magnetic field. In Section 4 we discuss the 
magnetic spectrum that is generated by IGW impinging on 
a magnetic field, while in Section 5 we discuss the wave en- 
ergy reaching the deep radiative interior. The effect of these 
waves on the development of the mean fiow is discussed in 
Section 6. We conclude with a discussion of the wave-field 
interaction. 



2 NUMERICAL MODEL 

We solve the coupled set of nonlinear magneto- 
hydrodynamic (MHD) equations in the anelastic ap- 
proximation. The anelastic approximation is valid when the 
fiow velocities and Alfven speed are much smaller than the 
sound speed. This approximation effectively filters sound 
waves, allowing a larger numerical time step. The equations 
are solved in a two-dimensional (2D) domain representing 
an equatorial slice of the Sun that spans the radius range 
from O.O5R0 to O.95R0. The details of the equations can 
be found in Rogers & Glatzmaier (2005) and Rogers & 
MacGregor (2010), which we will hereafter refer to as PI. 
The fundamental difference between this paper and PI, 
is that here we self-consistently solve for the convection, 
so that there is no need to artificially force waves on the 
top boundary. The radiative region extends from O.OSR© 
to O.TIR© and the convection zone extends from O.Tli?© 
to O.9OR0. For numerical purposes we impose another 
stable region on top of the convection zone, extending from 
O.9Oi^0 to O.93R0. The reference state thermodynamic 
variables are taken from the standard solar model by 
J.Christensen-Dalsgaard et al. (1991). 



The equations are solved in cylindrical coordinates. 
The radial dimension is discretized using a finite differ- 
ence scheme. The variables are expanded as a Fourier se- 
ries in the longitudinal dimension. The grid-space resolu- 
tion for all the models is 2048 longitudinal points and 1500 
radial points, with 1000 radial levels dedicated to the con- 
vection zone, convective-radiative interface and the magne- 
tized region. Solutions are time-advanced using the Adams- 
Bashforth method for the nonlinear terms and the Crank- 
Nicolson method for the linear terms. The velocity bound- 
ary conditions are impermeable and stress-free. The ther- 
mal boundary conditions are isothermal and the magnetic 
boundary conditions are perfectly insulating. The model is 
parallelized using Message Passing Interface (MPI). 

As described in PI, the magnetic field is written in 
terms of a vector potential, resulting in an evolution equa- 
tion for the vector potential. In all magnetic models we spec- 
ify a magnetic field between O.68R0 and O.7IR0, with the 
following form: 



A{r) = Aocxp 



(1) 



where A(r)z represents the vector potential, so that the 
magnetic field (B = V x A) is purely toroidal and is ini- 
tially confined between O.68R0 and O.7IR0, switching sign 
at O.695R0. Ao is varied substantially, leading to peak field 
strengths between 200 and 2 x 10^ G for the various models 
(see Table 1).^ 

The control parameters in these models are the mag- 
netic field strength, specified by Ao, and the thermal (k,), 
viscous (ly) and magnetic (77) diffusivities. The magnetic dif- 
fusivity is specified to be constant in radius, while both the 
thermal and viscous diffusivity are functions of radius with 
the following form: 



16crT(r) 



3p(r) xcp 



p{r) 



(2) 
(3) 



where cr is the Stefan-Boltzman constant, T(r) is the refer- 
ence state temperature, p{r) is the reference state density, x 
is the opacity and Cp is the specific heat at constant pressure. 
Therefore, we vary n and v by varying nrn and Vm- Because 
the thermal diffusivity fixes the heat fiux through the sys- 
tem and hence, the convective driving, for most models we 
keep Hirn (and hence, n) fixed. To investigate the effect of 
diffusivities and Prandtl numbers, we therefore vary Vm and 
77. For all magnetic models we run a non-magnetic model 
for comparison, a list of the models run and their control 
parameters (diffusivities) are presented in Table 1. 



3 THE INTERACTION OF INTERNAL 

GRAVITY WAVES AND MAGNETIC FIELD 

To carry out this study we start with a fairly evolved model 
from Rogers & Glatzmaier (2005), and impose the field de- 
scribed in Equation 1, as seen in Figure 1. The field is im- 

^ While the peak field strengths vary substantially, one should 
keep in mind that, in every model, the field strength varies from 
to the peak value in the form described in equation 1. 
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Table 1. Model parameters. Aq is the amplitude described in equation 1. Note that values of 1.66 x 10^^ correspond to 
peak magnetic field strengths of ^ 10^ G. u and 77 are the thermal, viscous and magnetic diffusivities listed in units 
of xlO^^cm-^s"^, at a radius of O.69R0, equivalent to mid-depth of the imposed magnetic field (where applicable). 
Models designated NF* refer to models without a magnetic field, models KV* similarly have no magnetic field but 
vary k,. Models designated BV* vary magnetic field strength keeping diffusivities constant. Models designated EV* 
vary 77, and models designated NV* vary u. 



Model 


Ao{10^^)iGauss/cm) 






r]{10^^cm^/s) 




T\.T A 


64 


23.6 


TVT A 
i\ A 


NF2 


NA 


64 


18.9 


NA 


NFS 


NA 


64 


14.2 


NA 


NF4 


NA 


64 


9.5 


NA 


NF5 


NA 


64 


4.7 


NA 


BVl 


1.660 


64 


23.6 


10. 


BV2 


0.166 


64 


23.6 


10. 


BV3 


0.0166 


64 


23.6 


10. 


BV4 


0.00166 


64 


23.6 


10. 


EVl 


1.66 


64 


23.6 


1. 


EV2 


1.66 


64 


23.6 


20. 


EV3 


0.166 


64 


23.6 


1. 


EV4 


0.166 


64 


23.6 


20. 


NVl 


1.66 


64 


18.9 


10. 


NV2 


1.66 


64 


14.2 


10. 


NV3 


1.66 


64 


9.5 


10. 


NV4 


1.66 


64 


4.7 


10. 


KVl 


NA 


90 


4.7 


NA 


KV2 


NA 


128 


4.7 


NA 



posed in an already evolved model because the initial time 
required to get convection started and waves set up is long 
enough that an initially imposed (two-dimensional) field will 
have decayed by some amount. The drawback of this method 
is that the field is placed in a region which is not quiescent 
and therefore immediately feels the presence of small scale 
structure. In these simulations it is impossible to disentan- 
gle "wave" from any other small scale motion. We therefore 
use the term "wave" to mean the following: if a variable is 
Fourier decomposed in longitude (as we do in these simula- 
tions), "wave" is any horizontal wavenumber, k, not equal to 
zero. Mean flow therefore, refers to horizontal wavenumber, 
k, equal to zero. 

3.1 Waves in the absence of a Magnetic Field 

One of the central questions with regard to IGW transport 
and mixing is the spectrum and amplitude of waves gen- 
erated at the convective-radiative interface. While there is 
still some debate, previous numerical calculations which in- 
clude the most relevant physics, but which are unfortunately 
limited to 2D (Rogers & Glatzmaier 2005), have predicted a 
broad spectrum of IGW generated at the base of the convec- 
tion zone. Indeed, the spectrum is very similar to that of con- 
vection. Deeper in the radiation zone the spectrum is dom- 
inated by high-frequency, low-wavenumber waves. Standing 
modes appear to be short-lived. 

Figure 2 shows the spectrum of horizontal wavenumber, 
k, versus frequency^ for Model NFS (no magnetic fleld) at 



three radii, O.71i?0, O.7Oi?0 and O.GQi?©. The color tables 
represent kinetic energy density. One can see that motion at 
the base of the convection zone has its highest amplitudes at 
lower wavenumbers and frequencies; there is little decrease 
in amplitude until wavenumber ^ 200 and frequencies of 
^ 100 - 150/iHz. 

The frequency and wavenumber dependence of the spec- 
trum of motions at the base of the convection zone can be 
seen more clearly in Figure 3, which shows the wave en- 
ergy density {p{vr + v^)) as a function of frequency (a) and 
wavenumber (b) at a radius of O.7OR0, for a variety of mod- 
els. Figure 3a shows a strong dependence of wave generation 
on frequency, with the lowest frequencies generated with the 
highest amplitudes and the highest frequencies nearly three 
orders of magnitude lower in amplitude, over a frequency 
range of 100/iHz. In Figure 4 we show fits to these curves 
for a variety of horizontal wavenumbers for model NFS; we 
find that the frequency dependence of the velocity ampli- 
tude can best be described as decaying like Cc;"^'^, although 
we find substantial variance, both with frequency and hor- 
izontal wavenumber, with exponents ranging from ^-2.2 to 
-3.2. If we compare this spectrum to that in the convection 
zone, we find very little difference (other than amplitude, 
see Figure 4 dashed line).^ There is some evidence for two 
slopes in the spectrum, with high frequencies better fit by 
a slightly steeper power law of cj"^'^ and lower frequencies 
fit by cj"^'^ which may indicate different generation mech- 
anisms. However, there is enough variance of the frequency 
dependence with horizontal wavenumber that it is difficult 



Note that for most of the spectra shown the frequency resolu- 
tion is 0.97/LtHz. However in 3.2 we study the effect of decaying 
field and in that case the frequency resolution is reduced. 



^ The similarity of the spectrum to that in the convection zone is 
not surprising since continuity of pressure at the interface implies 
continuity in velocity. 
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to draw any firm conclusions about generation mechanism. 
Since the energy density varies as approximately uj~'^'^ and 
the group velocity varies as approximately cj, we estimate 
that the wave energy flux varies as approximately uj~^'^ . 

Clearly, the amplitudes of the motions driven at the 
base of the convection zone depend on frequency. However, 
in Figure 3b we see that the generation of wave energy 
shows a very mild dependence of excitation on horizontal 
wavenumber with a maximum variation of one to two or- 
ders of magnitude (and generally less) over a horizontal 
wavenumber range of 100. At low wavenumbers, the ampli- 
tudes increase with horizontal wavenumber, similar to that 
predicted analytically in GLS and KTZ. However, from hori- 
zontal wavenumber 20-50, the amplitudes are relatively flat, 
and beyond that the amplitudes decrease. Motions in the 
convection zone show a similar, relatively flat dependence of 
energy on wavenumber, at least out to horizontal wavenum- 
ber 100-200. 

There has been significant discussion about the genera- 
tion of IGW at the base of the solar convection zone. In both 
the analytic works by GLS and KTZ the spectrum is calcu- 
lated by considering the stress imposed at the convective- 
radiative interface by convective eddies with smaller scales 
described by a Kolmogorov cascade. Another likely, and 
probably dominant, mechanism for generating IGW is that 
of overshooting plumes for which there is no analytic pre- 
diction of the spectra. In these simulations, we have both 
stresses at the interface due to convective eddies, as well as 
convective overshoot by plumes. The spectrum generated in 
these simulations with the energy density varying as uj~'^'^ 
and non-monotonic variation with wavenumber, is unlike ei- 
ther GLS or KTZ, which both predict energy density varying 
as uj~^k. 

There could be many sources for the discrepancy be- 
tween the spectrum produced by these simulations and those 
of analytic models. For example, both KTZ and GLS assume 
a Kolmogorov cascade which describes isotropic, inhomoge- 
neous turbulence, which is not likely representative of the 
base of the solar convection zone. The discrepancy could 
be due to the two-dimensional nature of these simulations 
compared to the one dimensional description of those ana- 
lytic models. Nonlinear interactions between waves included 
here but neglected in those models could give rise to dif- 
ferent correlation scales and times and therefore, a different 
spectrum. However, we think the most likely source for the 
discrepancy is the inclusion in these models of overshoot- 
ing plumes. Unfortunately, in this model we can not discern 
(and possibly it is not useful to do so) which mechanism 
(stresses and deformation of the interface versus overshoot- 
ing plumes) is responsible for the spectrum observed. It is 
plausible that the limited extent of a plume incursion could 
give rise to a broader wavenumber spectrum and that the 
short timescale of such an event leads to higher frequencies. 
However, our spectrum is not completely different than that 
obtained from analytic models and therefore, it appears that 
the most likely scenario is that both mechanisms are rele- 
vant. 

In addressing the spectrum of waves generated we can 
look to experimental results in limited circumstances. Ex- 
periments of IGW generated by an overshooting plume (An- 
song & Sutherland 2010), have predicted that the frequen- 
cies generated correspond to a narrow range (0.6-0.8) of 



the constant buoyancy frequency, N. This range is simi- 
larly found with sinusoidal disturbances dragged along the 
boundary Aguilar k, Sutherland (2006). While these results 
are interesting and likely relevant, it is hard to make the 
comparison with the models presented here, or the Sun for 
that matter, as the buoyancy frequency varies substantially 
with radius.^ 

In figure 2 we can also see rapid dissipation of the 
waves over the small radial range shown, we can particu- 
larly see low frequencies and higher wavenumbers dropping 
in amplitude with radius, consistent with simple radiative 
dissipation.^. We can investigate the dissipation of waves 
as a function of radius more readily by looking at particu- 
lar frequency and wavenumber combinations as a function 
of radius. In Figure 5a the dissipation of waves with ra- 
dius is shown for three different models with color repre- 
senting either varying thermal diffusivities (top four pan- 
els) or varying viscous diffusivities (bottom four panels) 
and linetypes indicating different frequencies, all with the 
same wavenumber (k=10,a,e,c,g) or different wavenumbers 
all with the same frequency (10/iHz,b,f,d,h). The left hand 
panels (a-d) are not adjusted for the variation in generation 
amplitude, discussed above, while the right hand panels (e- 
h) are adjusted. We can see two things immediately: (1) 
the wave dissipation is fairly independent of v and n and 
(2) the wave amplitudes in the deep radiative interior are 
fairly independent of frequency. We are careful here to dis- 
tinguish that the wave amplitudes in the deep interior are 
independent of frequency because the wave dissipation is 
dependent on frequency. This apparent contradiction arises 
because wave generation is strongly dependent on frequency, 
so that lower frequencies are generated at higher amplitude, 
but are simultaneously dissipated much more strongly, lead- 
ing to very little difference in wave amplitude in the deep 
radiative interior across a large range of frequencies^. In- 
deed, if we correct for the generation amplitude we can see 
a clear frequency dependence in Figure 5 (e,g), with higher 
frequencies damped less than lower frequencies, as expected 
from simple radiative dissipation. 

We can calculate the dissipation length by measuring 
that radius at which the wave amplitude decays by a factor 
of e, this is shown in Figure 6. In that figure the damping 
length is shown as a function of frequency and various line- 
types depict different horizontal wavenumbers. We immedi- 
ately see that the linear description of waves is qualitatively 
recovered; higher frequencies have longer damping lengths. 
However, we do not find the severe uj^ dependence predicted, 
but instead find something close to uj^ (uj^'^ is shown in 
the figure). This is qualitatively understood by consider- 
ing that the damping length is proportional to the group 
velocity divided by some damping rate, as outlined in Ku- 
mar et al. (1999). In the case of the Sun the ratio uofN 
is quite small and therefore the wave path is mostly hori- 

^ If one were to take this seriously, the range of 0.6-0.8 N, with N 
varying between zero and a few hundred /iHz, could give a broad 
frequency range as seen here. 

^ Clearly, the dissipation here is stronger than is likely in stellar 
interiors because of our necessarily large diffusion coefficients 
^ Of course this is not strictly true. As can be seen in the figure 
significantly higher frequency waves have higher amplitudes in 
the deep interior 
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zontal. Therefore, one should consider the horizontal group 
velocity duj/dkh, given approximately by N/m, where m 
is the vertical wavenumber. On the other hand, dissipation 
is likely more dominant in the vertical direction where the 
gradients are largest. Therefore, the dissipation rate is given 
by hikl. This leads to a damping length which is given by 
^ uj^ /{N'^kf^Hc). We should note here that our calculated 
damping lengths only recover the quantitative behaviour 
in damping length for a limited range of frequency and 
wavenumber combinations. At very low frequencies and high 
wavenumbers the vertical wavelength is relatively small and 
likely not well resolved throughout the domain (in the deep 
radiative interior we resolve vertical wavelengths of approx- 
imately 10^ cm, which if linear theory holds, would require 
u/k to be greater than ^ 10~^ for N ^ 100//Hz). On the 
other hand high frequencies and low wavenumbers have rel- 
atively large vertical wavelengths and the WKB approxima- 
tion is not well justified.^. 

Another potential problem with calculating the damp- 
ing length for a variety of frequencies and wavenumbers is 
that for many frequency /wavenumber combinations the mo- 
tion does not appear "wave-like" until fairly deep into the 
radiative interior. This is also apparent as two distinct slopes 
for many frequency /wavenumber combinations seen in Fig- 
ure 5. Nearer the convection zone the slope is steeper than 
deeper into the radiative interior. This is because close to 
the convection zone the waves are generated with rather 
high amplitudes and a simple estimate of the linearity of 
the waves given by Cp/u ^ 1 does not hold and linear the- 
ory is likely not valid. In this region mode- mode coupling is 
efficient and wave energy can be transferred to shorter scales 
where it is dissipated rapidly. Deeper into the radiative inte- 
rior, wave amplitudes are diminshed and one can clearly dis- 
cern wave like structure. Here, wave- wave coupling in scale 
is less efficient, leading to dissipation roughly in agreement 
with simple predictions. 

While the viscous diffusivity plays little role in the dis- 
sipation of waves with radius within the radiation zone, it 
does affect the amplitudes of the waves generated. Lower val- 
ues of the viscous diffusivity lead to less dissipation of the 
convective motions and therefore, slightly larger convective 
velocities. This in turn, leads to slightly higher amplitude 
motions below the convection zone at all frequencies and 
wavenumbers, as can be seen in Figure 7 a,b,c. In particu- 
lar, lower values of the viscous diffusion lead to higher ampli- 
tudes at higher frequencies. This is likely due to more dom- 
inant nonlinear terms (accompanying the larger Reynolds 
number) which lead to more efficient wave-wave coupling. 
This is seen in Figure 7d where, for a given frequency, lower 
values of the viscous diffusivity show shallower slopes in hor- 
izontal wavenumber, indicating better transfer between var- 
ious scales. Moreover, nonlinear interactions between waves 
could change the correlation times of those waves explaining 
the change in the frequency dependence of the spectrum. 



3.2 Waves in the presence of a Magnetic Field 

In MacGregor & Rogers (2010) it was shown that the 
strength of the magnetic field-IGW interaction depends on 
the difference between the Alfven frequency and the IGW 
frequency. If one considers a vertically propagating wave 
travelling from region 1 with field strength described by an 
Alfven velocity, uai, to an underlying region 2 with field 
strength described by ua2, then the reflection coefficient is 
given by (equation (27) that paper): 



R 



1 + c 



with 



7 2 2 

mi ^uj^ — k'^u\-^ 



1712 /CJ 



(4) 



(5) 



^ For frequencies of lOOyLiHz and horizontal wavenumber 5, the 
ratio of \r{dN'^ /dr)/N'^ is larger than 0.01 for most of the radia- 
tive interior 



Here m2 and mi represent the vertical wavenumbers in the 
two regions, uj the IGW frequency, and k is the horizontal 
wavenumber. The reflection coefficient thus depends on the 
ratio of the field strengths in region 1 and 2 (or, in the 
case presented here, on the derivative of the field strength 
with radius) and the IGW wave frequency. It was shown 
in MacGregor & Rogers (2010) that for low field strengths 
the reflection coefficient is small, but non-zero for a wide 
range of IGW frequencies. On the other hand, for high field 
strengths the reflection coefficient is large (one) but for a 
narrower range of IGW frequencies. This is generally true for 
both steep and shallow field strength gradients. Since these 
results imply that reflection could occur for a range of field 
strengths and IGW frequencies, we first turn our attention 
to how the IGW spectrum is altered in the presence of a 
magnetic field. 

When looking at the absolute spectrum of waves at 
some radius below the convection zone it is difficult to dis- 
cern the small differences in amplitudes between different 
models. Therefore, in the following we will often show the 
ratio of spectra for different models. We show the ratio of 
the kinetic energy spectrum of waves for Model BVl (strong 
magnetic field) to that of NFl (no field), in Figure 8a and 
the ratio of the kinetic energy spectrum of waves for Model 
BV2 (weak magnetic field) to that of Model NFl (no field) in 
Figure 8b at a radius of O.GQR©. These two models have all 
other parameters equivalent. In this figure regions which are 
blue indicate regions in frequency and horizontal wavenum- 
ber (f,k) space where the ratio is close to zero, meaning there 
is little or no energy in that region for the magnetic case, 
relative to the non-magnetic case. Regions which are red rep- 
resent ratios of ^ 1, implying that there is as much or more 
energy in that region of (f,k) space in the magnetic case as 
the non-magnetic case. What is immediately obvious in 8a 
are the dark ridges. In Figure 8b one can see dark ridges, but 
these features are both less severe and spread over a larger 
range in frequency and wavenumber. This is similar to what 
is seen in the analytic results of MacGregor & Rogers (2010). 
Namely, strong magnetic fields lead to strong refiection in a 
limited region of (f,k) space, while weak fields lead to weaker 
refiection in a larger region of (f,k) space. 

As described in Section 2 we impose a magnetic field 
which is purely toroidal. In the absence of a dynamo, the 
field cannot regenerate itself, and we would expect this 
toroidal field to decay in a diffusion time given by Tdiff = 
D'^/t], where D is the depth of the magnetic layer. In our sim- 



6 T.M. Rogers and K.B. MacGregor 



ulations the field is distorted by overshooting motions and 
IGW and therefore, the timescale over which diffusion oc- 
curs is much shorter. We show this effect in Figure 9, which 
displays the square of the mean toroidal magnetic field as a 
function of time in our simulations (solid lines) for Models 
EVl (black), BVl (blue) and EV2 (red) and that expected 
for pure diffusion, with the specified diffusion rate (dashed 
lines) . Clearly, our field is decaying much more rapidly than 
expected for pure diffusion. In fact, if one were to calculate 
the "effective magnetic diffusivity" , ?7e//, for each of these 
models one finds that Model EVl has ?7e// ~ 10?7ei;i, Model 
BVl has ?7e// ~ 3?7ei;i and Model EV2 has r]eff ~ 2?7e^;2• 
The dependence of the enhancement on 77 is reflective of the 
fact that higher initial values of 77 represent fields which are 
not as well tied to flows as their lower 77 counterparts and 
the field therefore, responds less to flow imposed distortion. 

Because filtering of the spectrum is a strong function 
of field strength, and because the field strength is a func- 
tion of time, we need to bear in mind that the wave filtering 
described above, could also be a strong function of time. 
We demonstrate this in Figure 10, where we show wave fil- 
tering as a function of time. In Figure 10 a,d,g,j we show 
the ratio of kinetic energies for the magnetic model com- 
pared to the non-magnetic model over a timescale of ~ lO^s 
and therefore, with a frequency resolution of 1//Hz, in 
(b,e,h,k) the timescale is 5 x lO^s and in (c,f,i,l) the timescale 
is 2.5 X lO^s, with the corresponding decrease in frequency 
resolution. Various rows show the ratio for different models 
(as outlined in the figure and figure caption). It is immedi- 
ately clear that the decaying field has a profound effect on 
the filtering of waves. In all cases, when the field is relatively 
in tact (early in the simulation) wave filtering is severe. For 
most models, filtering is reduced as the field decays. The 
changing properties of wave filtering are likely due to two 
effects: (1) the field amplitude decays and therefore the re- 
fiection coefficient described in equation 5 above decreases 
and (2) the wave-field interaction becomes less severe as the 
field is modified by the overlying convection. 

The field evolution for Model NV4 is shown in Figure 
11. The initial toroidal field does not remain stationary for 
long. Waves and convective overshoot impinge on the field 
causing smaller scale structure. Eventually the field is di- 
rectly impinged by overshooting plumes which entrain the 
field into the convection zone. Once there, the dynamics is 
best described as magnetoconvection. Field in the convec- 
tion zone is subsequently pumped back down and this very 
non-linear process continues, complicating the simple field- 
IGW interaction picture. Of course, once the field is dis- 
torted, the field gradient and thus the refiection coefficient 
become spatially and temporally dependent. A spatially de- 
pendent refiection coefficient could then allow transmission 
at one longitude while refiecting a wave at another. 

The above discussion describes models in which the 
magnetic Prandtl number, Prm (=z^/?7), is greater than one. 
In the solar tachocline Prm is less than one. To study the de- 
pendence of wave filtering on Pr^n we keep 77 fixed and lower 
V. The thermal and viscous diffusivities change substantially 
over the computational domain, as can be seen from their ra- 
dial forms in equations 2 and 3, so when we refer to Prandtl 
(magnetic and non-magnetic) numbers less than one we are 
referring to the region from O.GOR© and O.TIR©. All models 
have Prandtl numbers less than one, and we ran one model 



for which the magnetic Prandtl number is less than one. If 
we compare Figure 10 a,b,c with Figure 10 j,k,l we can see 
that as Prm is lowered below one, more wave energy is fil- 
tered in the field region. We can also see that Model BVl, 
shown in Figure 10 a,b,c is more susceptible to the effects of 
field diffusion on wave filtering when compared to 10 j,k,l, 
even though the magnetic diffusivity is equivalent in the two 
models, while the viscous diffusivity is substantially lower in 
NV4. This indicates the coupled nature of field and flow and 
their dependencies on all diffusivities. 

The decrease in kinetic energy in particular regions of 
(f,k) space indicates that incoming wave energy is reflected 
and ducted between field and no-field regions or between 
neighboring field regions (as implied in MacGregor & Rogers 
(2010)). In addition, it is likely that incoming IGW undergo 
mode conversion to Alfven waves; indeed, once in the field 
region gas motion is of mixed gravity- Alfven type. We now 
turn our attention to the magnetic spectrum that is gener- 
ated by the impinging IGW generated by convective over- 
shoot. 



4 MAGNETIC SPECTRUM 

Although it is possible that IGW are mode converted to 
Alfven waves, it is virtually impossible to distinguish that 
behaviour in these nonlinear calculations. Therefore, in the 
following we will discuss the spectrum of magnetic fiuctua- 
tions that develop from impinging IGW and how that mag- 
netic energy evolves in time, without specific regard to the 
exact manner in which IGW energy is converted to magnetic 
energy. 

Alfven waves are described by the simple dispersion re- 
lation UJA = VAk, where va is defined va = Bo/y/^irp. Be- 
cause the initial field we impose, Bo, varies substantially 
with radius, a range of Alfven speeds are possible. Fur- 
thermore, the convective motions drive a wide range of 
wavenumbers (see Figure 2), so one could expect a very 
broad range of Alfven frequencies to be generated. Indeed, if 
Alfven waves were directly excited by IGW with the corre- 
sponding wavenumber, one would expect Alfven frequencies 
in the range < uja < Bok/ ^Aiip. This would lead to fre- 
quencies ranging from zero to ^ 14/iHz for Model BV3 and 
from zero to ~ 1400/iHz for Model BVl (with the caveat of 
field decay discussed above). However, the interaction which 
converts IGW to magnetic energy is nonlinear and therefore 
depends on a variety of factors and direct excitation does 
not occur. 

As expected from the basic relationship above, we find 
that models with lower initially imposed field strengths gen- 
erally produce waves with lower frequencies and amplitudes, 
as seen in Figure 12. We see that large (r) amplitude mag- 
netic energy is concentrated at progressivley lower frequen- 
cies and wavenumbers for lower field strengths. Also note- 
worthy is broadband magnetic energy seen at all frequen- 
cies, due to nonlinear interactions. The models in Figure 
12 all have Pvm larger than one. Lowering Prm below one 
also affects both the kinetic energy spectrum and the mag- 
netic energy spectrum, as seen in Figure 13. Lower magnetic 
Prandtl number models show a high frequency branch not 
seen in higher Pr^ cases. When inspecting the ratio of the 
kinetic spectrum shown in Figure 13a we see that there is 
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a large region of (f,k) space for which the magnetic model 
does not support (kinetic) motion. This is the same region 
for which there appears to be decreased magnetic energy. 
This could indicate mode conversion, whereby when an in- 
coming IGW is reflected, the field responds to the reflection 
at a different wavenumber and frequency. However, this be- 
haviour is not completely understood. 

When wave motion is reflected at a given radius in 
the magnetized region, the now upward propagating wave 
is likely to be re-reflected at the top of the magnetized re- 
gion, when propagating from magnetized to unmagnetized 
layers, and so on. The energy in this ducted region could 
be dissipated, transferred to magnetic energy (as dicussed 
in this section), or transferred to mean flow as shown in 
Rogers & MacGregor (2010). Whatever the case, the energy 
propagating into the deeper radiative interior is affected. 



5 WAVE TRANSMISSION 

In the above discussion we showed that a toroidal magnetic 
field below the convection zone could filter the waves prop- 
agating in that region. In the above we investigated how 
this filtering affects the wave spectra, looking in detail at 
what frequencies and wavenumbers are filtered, how much 
they are filtered, and how that filtering depends on field 
strength and magnetic Prandtl number. Now we turn to the 
problem of how that filtering affects the amount of wave en- 
ergy reaching the deeper radiative interior. This is a critical 
consideration because in discussions of the ability of IGW to 
transport angular momentum or mix species, two main ques- 
tions arise: (1) the spectrum of waves generated and (2) the 
amplitudes of waves generated, i.e. how much wave energy 
is available. In the section 3.1 above and in Rogers & Glatz- 
maier (2005) we addressed (1) and (2) for non-magnetized 
models, in section 3.2 we addressed how a magnetic field 
alters (1) and here we will address (2) in the magnetized 
models. 

Figure 14 shows the integrated kinetic energy in waves 
from the bottom of the computational domain to O.GTR©, 
just beneath the position of the initially imposed magnetic 
field. What is shown in the figure is the ratio of integrated 
kinetic energy in waves in the magnetic case to the inte- 
grated kinetic energy in waves for the non-magnetic case. A 
value of one would indicate that the magnetic models pro- 
duce as much integrated wave energy in the deep interior as 
a non-magnetic case, while a value less (greater) than one 
indicates that the magnetic models produce less (more) en- 
ergy in that region. Initially (up to about 2 x 10^ s) the ratio 
for all models oscillates but stays very close to one. How- 
ever, after this initial period (which is approximately the 
time it would take many of the IGW to propagate through 
the magnetic layer) there is substantial variation with ratios 
reaching as low as 0.2 and as high as 4. 

Given the filtering seen in the above sections it is not 
surprising that the wave energy below the field would be re- 
duced from the non-field case. However, the fact that there 
is a substantial period during which the ratio is larger than 
one is unexpected. This can be seen more clearly in Figure 
15, which shows the kinetic energy in waves as a function 
of radius and time for the magnetic case (top) and the non- 
magnetic case (bottom). In this figure one can clearly see the 



hallmarks of IGW propagation with upward propagation of 
phase and downward propagation of energy. The top and 
bottom figures look virtually identical for the first 2 x 10^ s 
(which corresponds to the initial oscillatory behaviour seen 
in Figure 14), but then one can see a clear change in the field 
case with subsequent higher amplitude disturbances propa- 
gating inward in the magnetic model. One interpretation of 
this behaviour is tht the field, by virtue of reflection (and 
subsequent reflection at the base of the convection zone) 
traps wave energy in the field regio, causing the integrated 
energy in the deep radiative interior to drop below one. Once 
the field is sufficiently decayed or distorted the once trapped 
energy is allowed to propagate into the deeper radiative inte- 
rior, causing the integrated wave energy in the deep radiative 
interior to increase again, as seen in Figure 14. This process 
could be very time dependent given the spatial and tempo- 
ral dependene on the field in the region. This interpretation 
is somewhat supported by figure 14b, which shows the inte- 
grated wave energy ratio for model NV4/NF5 beneath the 
initially imposed field (as seen in figure 14a, solid line) as 
well as the ratio of the sum of the integrated wave kinetic 
and magnetic energy ratio in the field region. If the above 
interpretation is correct one would expect that the energy 
in the field region would be maximum when the energy in 
the deeper interior is minimum. This is true during much 
of the time shown in Figure 14b^, making this a plausible 
interpretation and one which is consistent with the theories 
laid out here and in (MacGregor & Rogers 2010; Rogers & 
MacGregor 2010). 

In Figure 15 one also notices a great deal of complexity. 
Even in the deep radiative interior, where the wave ampli- 
tudes have been severely attenuated, there is still consid- 
erable interaction leading to significant time and position 
variation of amplitudes and phase and group speeds. 



6 THE EFFECT OF THE MAGNETIC FIELD 
ON THE MEAN FLOW 

One of our main concerns when undertaking this study was 
how an imposed field might affect the ability of IGW to 
transport angular momentum in to the deep radiative inte- 
rior. In Figure 16 we show the ratio of the kinetic energy in 
the mean flow for the magnetized model compared to the 
non-magnetic model, in the field region (top) and in the 
deeper radiative interior (bottom), for various models. One 
can see that there is very little change in the mean flow in 
the deep radiative interior. This could be physical, in the 
sense that an imposed magnetic field does not affect the 
correlations which contribute to the mean fiow. However, 
it is possible that this is a numerical limitation. The mean 
flow in the deep radiative interior generally takes a substan- 
tial amount of time to grow. We are limited in how long we 
can run these models because of the field decay discussed in 
section 3.2; if we run these models much longer they would 
eventually revert to a non-magnetic case and the compari- 
son would be moot. So, while this is an important question 
to address, we are unable to do so definitively with these 
models. 



Although, it is clearly not true all of the time. 
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In the field region where the timescale for mean flow 
growth is likely to be much shorter (because the amplitudes 
of the waves are larger and because the region is shallower), 
we see that in most models the mean flow changes little from 
the unmagnetized model. However, in a few models we see 
that the mean flow grows and in our (arguably) most real- 
istic model (lowest Prm, largest field strength) we see that 
the mean flow grows substantially. This is similar to what 
was seen in (Rogers & MacGregor 2010)^. This is due to the 
fact that more wave energy is confined to the region, either 
by the ducting proposed in MacGregor & Rogers (2010) or 
by secondary excitation of IGW by field shown in Figure 15. 
More wave energy at a variety of wavenumbers and frequen- 
cies is likely to lead to more wave-wave interaction which 
could contribute to the mean flow. 



7 DISCUSSION 

In this paper we have continued our investigation of the in- 
teraction of IGW with magnetic fields. This paper addresses 
the complex problem of a wave spectrum self-consistently 
driven by an overlying convection zone interacting with an 
evolving magnetic field, similar what is expected in the solar 
interior. As expected from our previous work we find that 
magnetic fields reflect IGW in a way which depends on the 
magnetic field strength. This wave- field interaction excites 
magnetic fluctuations which also depend on field strength in 
a way which is consistent with the simple Alfven relation. 

Naively, the reflection of waves by a magnetic field 
would imply that less wave energy would make it to the ra- 
diative interior below the magnetic field. However, we find 
that this is not always the case. The wave energy below the 
field varies substantially in time and can be both larger and 
smaller than than in a non-magnetic case. 

Although one of our main motivations for this work was 
to understand how a magnetic field might alter the angular 
momentum transport in the deep radiative interior, we are 
cautious about conclusions drawn from these simulations for 
the mean flow in the deep radiative interior. This hesitation 
is because the timescale for the mean flow to evolve in the 
deep radiative interior is likely to be longer than the duration 
of these simulations. In the region of the field we see that, 
generally, the mean flow is very similar in the magnetized 
and un-magnetized models. However, in cases with lower 
values of vq we see an enhanced mean flow, due to increased 
wave energy (by reflection or secondary generation) in the 
region. 

Although many of our results agree with the analytic 
theory in MacGregor & Rogers (2010) and the simple mod- 
els in Rogers & MacGregor (2010), we find many complica- 
tions in these models. First, the field evolves substantially 
in time so that the simple picture of a stationary field with 
waves impinging on it is not an accurate description of the 
dynamics. The field evolution is dominated by the distortion 
of the field into smaller scales and the entrainment of field 
into the convection zone, and therefore is not diffusion dom- 
inated. This evolution affects wave filtering and, through its 



^ Although in that case it occurred even for Pr^^^ =1. 



affect on the convection, wave generation. These complica- 
tions (and more) are very likely to exist in the solar interior. 
The position of the toroidal field thought to give rise to the 
sunspots is in a region accessible to convective overshoot 
in order to allow for magnetic pumping. Therefore, distor- 
tion and entrainment of the field is likely. Additionally, the 
gas motions in that region are likely not linear waves, but 
some combination of nonlinearly interacting waves and con- 
vective overshoot as we noted in section 3.1. Second, even 
in the deep radiative interior, where wave amplitudes are 
low there appears to be substantial wave-wave interaction 
which leads to dissipation rates not well described by lin- 
ear theory. Finally, the secondary effect of field generated 
IGW (or, alternatively, time delayed and amplified waves) 
make estimates of wave energy in the deep radiative interior 
highly uncertain, given the unknown field strengths in the 
region and their evolution in time. Of course, all of this could 
be complicated even more by the presence of a primordial 
poloidal field, which will be the subject of a future work. 
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Figure 1. Time snapshot of the temperature perturbation, with 
black representing cold material and white reprsenting hot mate- 
rial, with the initial field depicted by the overlaid lines. Blue and 
yellow lines represent oppositely directed toroidal field. 
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Figure 2. Kinetic Energy Density spectrum for Model NFS at 3 
radii. The propagating waves seen in these models have lower fre- 
quencies than typical standing wave modes, which are short-lived 
and not seen here. Over a very short depth, high wavenumber/low 
frequency waves are efficiently dissipated. 
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Figure 3. Wave Generation. This figure shows the kinetic en- 
ergy density, calculated at O.TOR© as a function of frequency for 
wavenumber 10 (a) and wavenumber for a 10/iHz wave (b). Differ- 
ent colors represent different models: orange is model KV2, cyan 
is KVl, black is NFS, blue is NF4 and red is NF2. This figure 
shows the relative independence of wave generation on scale (b) 
and the strong dependence of that generation on frequency (a). 
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Figure 4. Wave Generation. Shown are wave kinetic energy den- 
sity as a function of frequency for various wavenumbers, k at a 
radius of 0.69 R© (solid line, in the Radiation zone) and a radius 
of 0.80 R0 (dashed line, in the Convection zone). The frequency 
dependence can be fit well by a power law with exponent around 
-2.7 to -3.0 (dotted line). Wave energy follows the distribution of 
convective energy quite well with very little difference except in 
amplitude. 
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Figure 5. Wave Dissipation. Left hand panels show wave kinetic 
energy as a function of radius, with amphtudes at the base of the 
convection zone as calculated, while the right hand panels show 
the wave kinetic energy as a function of radius, with amplitudes 
at the base of the convection zone adjusted to be equivalent, (a), 
(b), (e) and (f) show wave dissipation as a function of thermal 
diffusivity, ^c, with black lines representing model NFS, blue rep- 
resenting KVl and red representing KV2. Different linetypes in 
(a) represent different frequencies with solid lines represent 5/iHz, 
dotted representing 20)LiHz and dashed representing 80/iHz, all 
with horizontal wavenumber 10. Different linetypes in (b) rep- 
resent different wavenumbers, with solid lines representing k=3. 
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Figure 6. Wave dissipation length. Damping length measured 
in the simulations as the radius at which the wave amplitude 
dropped by a factor of e, for wavenumbers 15 (solid), 18 (dot- 
ted), 22 (dashed), 23 (dash-dot), 25 (dash-triple dot) and 28 (long 
dash). The straight line represents a power law of cj^-^. 
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Figure 7. Kinetic energy density spectrum for models NFl (a), 
NF3 (b) and NFS (c) at CTOR©, showing that lower values of 
u lead to the generation of higher frequency waves, (d) shows 
the amplitude of kinetic energy density as a function of horizon- 
tal wavenumber, for f=50/iHz, with the red line indicating model 
NFS, green line representing model NF3 and the blue line repre- 
senting model NFl. 
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Figure 8. The ratio of kinetic energy density for models 
BVl/NFl (a) and for models BV2/NF1 in (b) at a radius of 
O.69i?0. Dark blue represents a ratio of zero, indicating there is 
substantially less kinetic energy in that region of (f,k) space in 
the magnetic case than in the non-magnetic case, while red rep- 
resents a ratio of 1.2 indicating more energy in that region of 
(f,k) space in the magnetic model compared to the non-magnetic 
model. This clearly demonstrates field strength dependent wave 
filtering. 
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Figure 9. The energy in the mean field as a function of time 
for models with varying values of the magnetic diffusivity 77, EVl 
(black), BVl (blue) and EV2 (red). Dashed lines represent the 
diffusion expected for dissipation of the large scale field calcu- 
lated assuming TfUff = /rj^ where D is the initial depth of the 
magnetic layer. 
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Figure 10. Wave filtering in the presence of field. This image 
shows kinetic energy density ratio spectra for different models 
over different times. Different rows represent filtering for different 
models as indicated in the figure. Different columns represent 
the filtered spectrum over different time intervals with the left- 
most plots (a,d,g,j) representing integrations over lO^s, (b,e,h,k) 
representing integration over 5 x lO^s and (c,f,i,l) representing 
integrations over only 2.5 X lO^s. One can see that at early times, 
when the field has not decayed, the filtering is strong in all cases. 
However, after the field has decayed more (left most plots) the 
effect of filtering has decreased. Note that all of these models 
have the same (strong) initial field strengths, but varying degrees 
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Figure 11. Time snapshots of the magnetic field evolution. The 
top panel shows the toroidal field, while the bottom panel shows 
the radial field, both progressing in time to the right. The first 
snapshot of the radial field shows no structure as the field ini- 
tially imposed is purely toroidal. Radial field is created in time 
as overshooting plumes distort and entrain the field. 
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Figure 12. Magnetic Energy spectrum at O.TOR© as a function 
of field strength, with (a) showing the lowest field strength and 
(c) showing the highest field strength. As expected from the sim- 
ple dispersion relation for Alfven waves the highest field strength 
supports the highest frequency waves. 



22 r.M. Rogers and K.B. MacGregor 




Figure 13. Wave filtering in the presence of field. On the left we 
show the ratio of kinetic energy in the field case to the non-field 
case (same as what is shown in figure 8), on the right we show the 
magnetic energy spectrum. The circled areas show regions where 
there is little kinetic energy in the magnetic cases and similarly, 
on the right, where there is little magnetic energy in the magnetic 
spectrum. In this model the magnetic Prandtl number is less than 
one and the field is sufficiently tied to the ffow so that the region 
of (f,k) space where there is little kinetic energy is also that region 
of (f,k) space where there is little magnetic energy. 
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Figure 14. Integrated wave energy in the deep radiative interior. 
In this figure we show the ratio of kinetic energy density for the 
magnetic case to the non-magnetic case integrated from the base 
of the computational domain to O.67R0. Different hnetypes repre- 
sent various models: solid line NV4/NF5, dashed line NV3/NF4, 
dotted line NV2/NF3, dash-dotted line NV1/NF2, dash-three dot 
line BV1/NF2 and long dashed line BV2/NF2. Initially the mod- 
els all oscillate around one, but over time the models diverge with 
substantial excursions ranging from substantially less energy in 
the interior in the field case to substantially more.(b) The same 
ratio shown in (a) NV4/NF5 (solid line) and the sum of kinetic 
and magnetic wave energy ratio for NV4/NF5 in the field region, 
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Figure 15. IGW generation by magnetic field. This figure shows 
the integrated wave energy as a function of time (x-axis) and ra- 
dius (y-axis) beneath the magnetic field. On the bottom we show 
the unmagnetized Model NF5, while on the top we show the mag- 
netic Model NV4. Initially (to about 2xl0^s) the two models are 
identical, subsequently one can see that the magnetic model pro- 
duces more kinetic energy in the deep interior. The circled region 
shows increased kinetic energy in a group of IGW. Subsequent 
motion in the magnetized model also has larger kinetic energy. 
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Figure 16. Integrated mean flow in the radiative interior. In this 
flgure we show the ratio of the kinetic energy in mean flow for 
the magnetized case to the unmagnetized case for various models, 
sohd line represents NV4/NF5, dotted line is NV3/NF4, dashed 
line is NV2/NF3, dash- dot Une is BVl/NFl and dash-three dot 
line is EVl/NFl. In the top panel we show this ratio in the field 
region (integrated from O.GTR© to O.TOR©), and in the bottom 
panel we show this ratio below the field (integrated from the 
bottom of the domain to O.GTR©). Below the field we see that this 
ratio is nearly one, with some minor excursions above and below 
one. However, in the field region we see that both the model with 
the lowest u (NV4/NF5) and that with the lowest r] (EVl/NFl) 



